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This paper presents a new technique for achieving spectral accuracy and fast computational 
performance in a class of black hole perturbation and gravitational self-force calculations involving 
extreme mass ratios and generic orbits. Called spectral source integration (SSI), this method should 
see widespread future use in problems that entail (i) point-particle description of the small compact 
object, (ii) frequency domain decomposition, and (iii) use of the background eccentric geodesic 
motion. Frequency domain approaches are widely used in both perturbation theory flux-balance 
calculations and in local gravitational self-force calculations. Recent self-force calculations in Lorenz 
gauge, using the frequency domain and method of extended homogeneous solutions, have been able 
to accurately reach eccentricities as high as e ~ 0.7. We show here SSI successfully applied to Lorenz 
gauge. In a double precision Lorenz gauge code, SSI enhances the accuracy of results and makes a 
factor of three improvement in the overall speed. The primary initial application of SSI for us its 
raison d’etre-is in an arbitrary precision Mathematica code that computes perturbations of eccentric 
orbits in the Regge-Wheeler gauge to extraordinarily high accuracy (e.g., 200 decimal places). These 
high accuracy eccentric orbit calculations would not be possible without the exponential convergence 
of SSL We believe the method will extend to work for inspirals on Kerr, and will be the subject of a 
later publication. SSI borrows concepts from discrete-time signal processing and is used to calculate 
the mode normalization coefficients in perturbation theory via sums over modest numbers of points 
around an orbit. A variant of the idea is used to obtain spectral accuracy in solution of the geodesic 
orbital motion. 
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I. INTRODUCTION 

Merging compact binaries are a promising source of de¬ 
tectable gravitational waves. Accurate theoretical mod¬ 
els serve as templates to assist detection and will aid in 
estimating an event’s physical parameters. Three com¬ 
plementary theoretical approaches exist 1] for comput¬ 
ing relativistic binaries: numerical relativity UlSj, post- 
Newtonian (PN) theory [U 5], and gravitational self¬ 
force (GSF) and black hole perturbation (BHP) calcula¬ 
tions 03 (Ml- The effective-one-body (EOB) formalism, 
drawing calibration of its parameters from the above ap¬ 
proaches, then provides a synthesis HDHT51- 

The GSF approach assumes the existence of, and ex¬ 
ploits, a small ratio q = p/M <C 1 between the compo¬ 
nent masses. The field and motion of the smaller body 
are calculated in a perturbation expansion in powers of 
q mm- Though restricted to small q, the GSF is valid 
throughout the strong field regime. GSF/BHP calcu¬ 
lations are most relevant to potential future eLISA ob¬ 
servations of extreme-mass-ratio inspirals (EMRIs) q ~ 
10~ 7 -10~ 4 T8\ but might pertain to Advanced LIGO 
observations if there exists a fortuitous population of 
intermediate-mass-ratio inspirals (IMRIs) q ~ 10~ 3 -10~ 2 
tin m- The dominant approach to the GSF treats 
the small body as a point mass [8-, then calculates the 
metric perturbation and the local self-force using mode- 
sum regularization m- Calculations are done directly 
in the time domain (TD) [22H25] or via decomposition 
into Fourier-harmonic modes in the frequency domain 
(FD) p?7H - [TTrjl . Alternative means of calculating the GSF 


include effective source calculations j3lk33| and direct 
Green function calculations [34H36] . 

The PN approach has no restriction on q but is most 
accurate for wide, low frequency orbits. Just as the GSF, 
PN, and NR approaches separately inform EOB, there 
has been considerable activity in recent years in making 
comparisons between GSF/BHP and PN theory J2B1 [3Tb 
[35) , including calculations at very high accuracies HOF 
l43j . These high precision calculations, until recently all 
done for circular orbits, utilize the analytic function ex¬ 
pansion formalism of Mano, Suzuki, and Takasugi (MST) 
i4J and make use of arbitrary precision coding. 

Several of us were intrigued by the idea of extending 
MST calculations, and these comparisons, to include ec¬ 
centric orbits. Initial results of that now successful ef¬ 
fort will be described elsewhere [25] , but the project led 
to the necessary development of the technique reported 
here. Modeling EMRIs with large eccentricities is essen¬ 
tial, since astrophysical considerations suggest m (203 
they have a distribution peaked about e ~ 0.7 [2j3j as 
they enter eLISA’s passband. Advanced LIGO inspirals 
are expected to have nearly circular orbits, but small ec¬ 
centricity corrections may be important m- 

The spectral technique described here benefits FD cal¬ 
culations of the “geodesic self-force” (i.e., first-order per¬ 
turbations derived using geodesics of the background ge¬ 
ometry) in E/IMRIs with eccentric orbits. In eccentric- 
orbit FD calculations the Fourier transform spreads the 
influence of the point particle source across a range of 
radii. Mode by mode, the resulting source functions are 
integrated against a Green function over this radial li- 
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bration region, a procedure that has been followed for 
decades [45 for BHP problems. FD calculations of ec¬ 
centric orbit GSF became feasible after Barack, Ori, and 
Sago jJS] found the method of extended homogeneous 
solutions (EHS), thus allowing Fourier synthesis at the 
particle location without encountering Gibbs behavior. 
Originally demonstrated for scalar models [12 and early- 
on extended to master equations in the Regge-Wheeler- 
Zerilli (RWZ) formalism [27] . EHS has recently been ap¬ 
plied to coupled systems in Lorenz gauge [29l 30]. 

To those familiar with EHS, the new method can be 
outlined briefly here. See Secs. m and m for details. 
Here we couch the discussion in terms of the RWZ case 


(Sec. Ill B), where EHS entails calculating arbitrarily 
normalized causal homogeneous solutions 
and integrating them in product with stress-tensor pro¬ 
jections over the source region. The result is a set of nor¬ 
malization coefficients Cff ln that encode the orbital mo¬ 
tion’s imprint in the field perturbation. Then extended 
homogeneous solutions are assembled in the TD, and sub¬ 
sequently abutted at the instantaneous particle location. 

The new technique provides a means of calculating the 
(or their equivalent in other gauges) with spectral 
accuracy. The integral for the normalization coefficients 
is typically manipulated [27, '49] [50] into a form like 


r ± — 
w Imn 
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where Ef^^t) is a periodic function of the radial mo¬ 
tion derived from spherical harmonic projection of the 
point source and integration over the Xj^ mn (r). Details 
are found in Sec. Ill but for the nonce it is enough to say 


that computing (1.1) is difficult to do with ODE or nu¬ 


merical quadrature integrators at high accuracies beyond 
double precision and is impossible to do at extraordinar¬ 
ily high accuracies like 100 or more decimal places. The 
new method, called spectral source integration (SSI), re¬ 
places the integral with a remarkably simple sum 


1 

Cf 1 =- (t k ) e inQrtk , (1.2) 

lmn NW lmn ^ lmnK k> 

which involves merely sampling the source function 
£ ; ^ n (f) at a modest number N of equally-spaced points 
around the closed radial motion. This sum converges 
exponentially with increases in N. 

The FD approach with use of Fourier series (FS) has 
been a part of BHP theory for decades. The FS and 
normalization coefficients converge exponentially with n, 
allowing the FS to be truncated. The new method makes 
a crucial use of that standard approximation, recogniz¬ 
ing that truncation of the FS representation (of e.g., a 
source term) generates a bandlimited function. That in 
turn invokes the machinery of the Nyquist-Shannon sam¬ 
pling theorem. The truncated FS can itself be replaced 
by discrete equally-spaced sampling of the TD function. 
Then, discrete sampling and periodicity allow a discrete 


function of finite length N to serve as an accurate TD 
representation. Furthermore, the finite discrete function 
is dual to a discrete Fourier transform (DFT) spectrum, 
computable with an FFT. The DFT spectrum is an ap¬ 
proximation, between its Nyquist frequencies, of the orig¬ 
inal FS spectrum, but can be made exponentially accu¬ 
rate with increases in N. It is then possible to replace 
integrals like (1.1) with finite sums like (1.2) and achieve 
spectral convergence there too. In essence, SSI provides 
a completion of the FD approach by bringing to bear 
concepts in discrete-time signal processing. 

This paper shows application of SSI to FD BHP 
and geodesic GSF calculations of eccentric Schwarzschild 
E/IMRIs in both RWZ and Lorenz gauges. We also 
demonstrate in Sec. [TT] that a related approach pro¬ 
vides arbitrarily accurate solutions of the geodesic equa¬ 
tions themselves. SSI may be applicable to Kerr BHP 
pun, 52: and GSF calculations, the subject of an upcom¬ 
ing paper. In addition, SSI has the potential to benefit 
the Green function approach to GSF calculations [36] . 

This paper is organized as follows. Sec. |H] considers 
the orbital problem. In Sec. IIA we review bound ec¬ 
centric geodesic motion about a Schwarzschild black hole 
and set the notation. Sec. |II C| describes the spectral ap¬ 
proach for integrating the orbit equations with geometric 
convergence, and shows numerical results. Appendix [A] 
gives a simple analytic calculation of the exponential fall- 
off in Fourier coefficients in part of the orbital problem. 
Next the SSI method is described in Sec. |III| through its 
application in the RWZ formalism to provide spectral 
solution of master equations. A brief review of how the 
RWZ problem is solved in the FD using EHS is given in 
Sec. IIIA Then Sec. III B| lays out the SSI method, the 
heart of this paper, and displays a set of numerical re¬ 
sults. We discuss some related findings in the numerical 
analysis literature in Sec. |III C| Having shown SSI ap¬ 
plied to a single perturbation equation, we present next 
in Sec. |IV| its application in Lorenz gauge, demonstrating 
that the method allows systems of equations to be solved 
with spectral convergence. Our conclusions are drawn in 
Sec. El 

In this paper we set G = c = 1 and use the metric 
signature +2. 


II. SPECTRAL INTEGRATION OF BOUND 
ORBITAL MOTION 

The new method is first applied to solving the equa¬ 
tions of bound geodesic motion. This proves to be a 
necessary first step to using SSI to solve the first-order 
perturbation equations when working at accuracies well 
beyond double precision. At double precision, it leads to 
a more efficient computation of the orbit. We consider 
geodesic motion about a Schwarzschild black hole in this 
paper. Application of SSI to general orbits about a Kerr 
black hole will be taken up in a subsequent paper. We 
begin with a brief review of the problem and notation. 
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A. Geodesic motion and the relativistic anomaly 

We consider generic bound motion of a small mass p, 
taken to be a point particle, around a Schwarzschild black 
hole of mass M in the test body (geodesic) limit g/M —> 
0. Schwarzschild coordinates x^ = ( t , r, 9, p) are used, 
with the line element having the form 

ds 2 = - fdt 2 + f~ l dr 2 + r 2 (d9 2 + sin 2 6 dtp 2 ) , (2.1) 

where f(r) = 1 — 2 M/r. 

Let the worldline of the particle be described by the 
functions z“(t) = \t p {r), r p (r), 9 p (t), P p (t)} of proper 
time r (or some other convenient curve parameter). Sub¬ 
script p indicates location of the particle. The four- 
velocity is u a = dXp/dr. Without loss of generality the 
motion is confined to the equatorial plane, 9 p (r) = ir/2. 

The orbit is parametrized in terms of the (dimen¬ 
sionless) semi-latus rectum p and the eccentricity e (see 
[48l 853]). These constants are related to the usual con¬ 
stant specific energy £ = —u t and specific angular mo¬ 
mentum C = u v . Additionally, pericentric r m ; n and 
apocentric r max radii are introduced, which are related 
to p and e by the following equations 


^max^min 

T max • min 

(2.2) 

M(r max + r m in )’ 

^max ^min 

pM 

^max — i 

1 — e 

pM 

^min — 

1 + e 

(2.3) 


Bound eccentric orbits satisfy £ < 1 and C > 2\/3 M. 
These in turn imply p > 6 + 2e, with the boundary of 
stable orbits p = 6 + 2e being the separatrix [48] . 

As is usual, r is replaced as the curve parameter by 
Darwin’s relativistic anomaly y, in terms of which the 
radial position is given a Keplerian-appearing form [54] 


To solve (2.51 and (2.61, each equation can be regarded 
as either a numerical quadrature or an initial value prob¬ 
lem (IVP) [SHI- Cutler, Kennefick, and Poisson [4S] took 
the former approach and used Romberg’s method. In 
the more complicated Kerr geodesic problem, Drasco and 
Hughes m initially solved for the motion using a nu¬ 
merical quadrature routine but later switched to use of a 
quasi-analytic approach developed by Fujita and Hikida 
[571 . (Indeed, this quasi-analytic method involving rapid 
evaluation of elliptic integrals stands as a third rou te to 
solution.) In more recent work [371 [211 [301 [S3] > Eqns. (2.5) 


and (2.61 have simply been integrated using Runge-Kutta 
routines. At double precision the distinction is trivial and 
errors in the orbit are of minimal concern. Recently, how¬ 
ever, several of us have turned attention [45] to making 
extraordinarily high precision (e.g., 200 decimal place) 
BHP and GSF calculations for eccentric EMRIs using 
the MST formalism [44] (henceforth the MST code). It 
proved necessary to develop a new means of efficiently 
calculating the orbit to arbitrary precision, as well as 
doing the same for the perturbation source integration 
(Secs. |HIB] and |TVB |. 

The MST code is written in Mathematica to make use 
of its arbitrary precision functionality. Initially, we used 
its NDSolve function to compute orbits but found such 
integrations became prohibitively expensive for errors of 
order < 10” 40 . The alternative approach we found turns 
out to be a simple application of the SSI concepts. In 
fact, the arguments laid out in the next two subsections 
are key to understandi ng th e SSI development. Shortly, 
we will discuss solving ( |2.5[ ) to obtain t p (x) (integration 
of (2.6) follows in like fashion). But first we address some 
general considerations. 


B. Spectral integration: general considerations 


r P (x) = 


pM 


(2.4) 


1 + e cos y 

The equations for the remaining functions take the form 


dt p = r p (x) 2 / ( P - 2) 2 - 4e 2 

d\ M(p— 2 — 2e cos x) V P ~ 6 — 2e cos y ’ 

dr p Mp 3 / 2 I p — 3 — e 2 

dy (1 + ecosy) 2 V p — 6 — 2ecosy ’ 

d<p P _ j V 
dx \ p- 6-2e cos y 


(2.5) 

( 2 . 6 ) 
(2.7) 


The last equation, describing azimuthal motion, has an 
analytic solution 


Ap(X') = \j 


4 p 

p — 6 — 2e 



4e \ 
~ p — 6 — 2e) ’ 


( 2 . 8 ) 


where F(x\m ) is the incomplete elliptic integral of the 
first kind 155] . The other two equations are typically 
solved numerically. 


Let dl/dx = g(x) with g( y) (the source) being both a 
periodic and a smooth function. We are interested in in¬ 
tegrating g to find J(y). We can assume g( y) is complex, 
but in orbital motion applications the functions will be 
real. The periodicity of g suggests utilizing a FS expan¬ 
sion and then calculating the integral for J(y) term by 
term. At first glance this approach is not very helpful 
since, even if we truncate the FS, the expression for J(y) 
would require computing a large number of definite inte¬ 
grals numerically for the FS coefficients. Fortunately, the 
smoothness of <?(x) helps in several ways. In many cases, 
the FS amplitudes Q n will fall in magnitude exponentially 
(shown numerically for orbital motion presently; see also 
Appendix A). Even in calculations with hundreds of dec¬ 
imal places of accuracy, the FS can then be truncated to 
a modest number of terms. At whatever adopted level of 
accuracy, replacing g(y) with a truncated FS introduces 
an approximation that is bandlimited. 

We then recall that bandlimited signals play a key role 
in the Nyquist-Shannon sampling theorem: a function 
that contains only frequencies / with |/| < B is com- 
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pletely determined by its discrete (equally-spaced) sam¬ 
ples (in this case in y) occurring at the Nyquist rate 2 B 
(i.e., with spacing Ay = |i? _1 ). If we combine discrete 
sampling with the periodicity of radial motion, then only 
a finite total number TV of samples in y need be con¬ 
sidered. We replace g( y) again-this time with its finite 
sampling g k = g(xk) = g(k Ay), where k = 0,...,TV - 1. 
This new representation of the source has its own DFT 
spectrum Q n (with n = —TV/2,... ,7V/2 — 1), which can 
be computed with the FFT algorithm !56j. In contrast to 
the FS spectrum Gm the DFT spectrum Q n exhibits a pe¬ 
riodicity of its own, Gn+jN = Gn, for arbitrary integer j. 
However, aliasing can be avoided if the DFT spectrum is 
only used at the TV frequencies within its Nyquist bounds. 
Then for an accuracy goal that is sufficiently high (i.e., 
high enough TV, found iteratively), the DFT spectrum 
G n is virtually indistinguishable from the FS spectrum 
G n - Using the DFT representation, it is then possible 
to compute g{ y) at any location either via Fourier in¬ 
terpolation or using the Whittaker cardinal function |55j 
on the circle (i.e., convolution with the Dirichlet kernel). 
Furthermore, the source can be integrated or differenti¬ 
ated term by term to accuracies comparable to the initial 
goal. 

To summarize: 

• The (perhaps complex) function g(y) is periodic 
and C°°. 

• It can be represented as a FS with spectrum Gn 
with n —> ±oo. 

• The FS spectrum can be truncated to some n m i„ < 
n < n max subject to an accuracy goal. 

• The approximate but very accurate truncated FS 
is a bandlimited function. 

• The Nyquist-Shannon sampling theorem implies 
the truncated FS representation can itself be re¬ 
placed in the TD with discrete sampling. 

• Sampling plus periodicity implies a discrete repre¬ 
sentation of finite length TV. 

• Finite sampling representation in the TD implies 
one-to-one correspondence via the DFT with a FD 
periodic spectrum Gn- 

• The DFT spectrum within the Nyquist range ap¬ 
proximates well the original FS spectrum if TV is 
sufficiently large, allowing Gn —► Gn- 

• The DFT representation in the TD can be inte¬ 
grated and interpolated to spectral accuracy. 


C. Spectral solution of the orbital motion 


In practice, the orbit equations (2.5) and (2.6) have 
source functions that are real and even. Hence we can 



FIG. 1. Equally spaced in y sampling of p = 50, e = 
0.7 orbit. The complete orbit is split into TV = 42 samples 
(Ay = 0.1496) and spectral integration requires only A f = 
22 points between y = 0 and y = 7r (inclusive) to achieve 
double precision accuracy. The values of dt p /d\ need only be 
calculated at the indicated points to provide double precision 
integration and interpolation anywhere on the orbit. 


represent them with a discrete cosine transform (DCT) 
m- In turn the integral for t p ( y) (for example) will 
be represented by a discrete sine transform (DST), with 
an additional term linear in y. Furthermore, the orbital 
source functions are not only periodic but have reflection 
symmetries across both periapsis (y = 0) and apapsis 
(y = 7 r). These symmetries narrow the form that the 
DCT can take to be either type I or II [8j0J. We utilize 
the type I (referred to as DCT-I) algorithm with unitary 
normalization (making the DCT-I its own inverse). 

In the general discussion above, we imagined dividing 
the entire orbit into TV intervals with Ay = 2n/N. For 
the DCT-I, this spacing is maintained and (assuming TV 
is even) the half orbit from y = 0 to y = 7r is split into 
TV/2 intervals. The DCT-I utilizes A f = TV/2 + 1 sample 
points by including the end points at both y = 0 and 
y = 7r. In terms of the number of samples the domain is 
split into J\f — 1 intervals. The locations of the samples 
are 


Xk ~JT^l' fee 0,1,...,AA-1. (2.9) 

Then at the A f points we denote the samples of the source 
function as g k = g{xk)- The (real) Fourier coefficients are 
given by 


Gn 




1 

2 


9o + 2 ( — 


( 2 . 10 ) 


M-2 

+ ^2 9k cos (rcyjfe) 
fc=l 
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FIG. 2. Number of sample points A f between \ = 0 and 
X = tt needed to represent dtp/dx = g(x) 1° a prescribed 
accuracy. The ratio of magnitudes of the smallest to largest 
Fourier coefficients of g(x) gives an estimate of the relative 
accuracy. The linear scaling of Af versus digits of accuracy 
indicates geometric fall-off in the spectral components of g{x). 
Away from the separatrix this relation is largely independent 
of p. 


FIG. 3. Number of sample points A f between \ = 0 and 
X = 7r needed to represent dt p /dx = g{x) as a function of ec¬ 
centricity e at fixed accuracy of 150 decimal places. Of course 
as the eccentricity approaches unity and the radial period be¬ 
comes infinite, so too does the number of needed samples. 
Still, for any astrophysically relevant orbit, we can sample 
the orbit to impressive accuracy with a modest number of 
points. 


Like the more general DFT, this expression is exactly 
invertible and we can recover the original samples in the 
X-domain 


Then, with T r in hand, the fundamental frequencies can 
be computed 


9k = 


Af — 1 




( 2 . 11 ) 


M -2 

+ ^2 Gn COS ( nxk ) 


We can then use the spectral amplitudes to provide a 
Fourier interpolation to arbitrary \ 


g(x) = 


AT- 1 


^Go + \Gu-i cos [(Af - l)x] 

AT -2 

+ Gn cos (n\) 

n= 1 


. ( 2 . 12 ) 


Integrating Eqn. (2.12) yields the sine expansion for the 
time 


t P (x) = 


AT-1 


, 1 n sin [(AT — l)x] 
2^oX+2^-i {Af _ 1) 


JT-2 


V —Gn sin (n\) 

z ' n 


n—1 


. (2.13) 


Having found t p (x) we can obtain the radial period T r 
from the leading Fourier amplitude Go 


T r = 


AT-1 


ttGo- 


(2.14) 




2tt 

T~’ 

1 - r 


&np — 


M 2?r ) 

T 
± r 


(2.15) 


From a practical perspective, the DCT-I can be com¬ 
puted with 0(Af\n.Af) speed using either the Fouri- 
erDCT function in Mathematica or the FFTW routine 
in C coding. Fig. [T] provides a picture of how efficient 
this method is. For this orbit we need only Af = 22 
samples to achieve double precision accuracy in the or¬ 
bit integration. In fact, all we need know are the source 
functions at the indicated points and we can interpolate 
to double precision accuracy anywhere in between. From 
a practical standpoint, we guess a value of Af and esti¬ 
mate the error by computing the ratio of the smallest to 
largest Fourier coefficients \Gm-i/Go\- If that ratio fails 
to meet our prescribed accuracy goal, we simply increase 
Af and repeat the procedure. Given that the DCT is so 
fast to compute, we are able to solve the orbit equations 
to hundreds of digits of accuracy within a few seconds. 
Fig- ID shows how remarkable and modest the scaling is 
in the number of required sample points Af as a function 
of prescribed accuracy. Fig. [3] shows how the number of 
needed sample points grows with increasing eccentricity 
(given a fixed accuracy goal). Even at very high eccen¬ 
tricities, e ~ 0.9, the number of samples is quite rea¬ 
sonable. Thus, with this approach the integration of the 
orbit becomes a trivial cost, even for astrophysically in¬ 
teresting eccentricities (e ~ 0.7) and high accuracy (MST 
code) applications. 
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III. SPECTRAL SOURCE INTEGRATION IN 
THE RWZ FORMALISM 


One of the principal goals of this paper is to describe 
our new means of applying spectral techniques (i.e., SSI) 
to integrate the source problem in black hole perturba¬ 
tion theory to high accuracy. In this section we show the 
simplest application of SSI, namely solution of master 
equations in the Regge-Wheeler-Zerilli (RWZ) formalism 
for generic orbits. Before detailing the SSI technique in 
Sec. |III B[ we first briefly review the now standard way 
m of solving master equations using FD decomposi¬ 
tion and the method of extended homogeneous solutions 
(EHS) HSj], and in the process set the notation. 


A. The RWZ formalism and EHS method 


We begin with a RWZ master equation in the TD 

(-^2 + ^2 - ^O - )) = S lm (t,r), (3.1) 

where r* = r + 2Mln(r/2M — 1) is the usual tortoise 
coordinate. Here Vj(r) is either the Zerilli potential (l + 
m even) or the Regge-Wheeler potential (l + m odd). 
The source contains terms proportional to the Dirac delta 
function and its first derivative 


the TD master functions over one radial period 


Xlmn (0 — 


T, 


dt T im (t,r) e 1 ' 


r J 0 


Zimn{r) = — / dt Sim{t,r)e l ‘ 


(3.6) 

(3.7) 


L r JO 


The master equation then takes on the following FD form 
d 2 


drl 


+ W 2 - y,(r) x lmn (r ) = Z lmn (r). (3.8) 


(Throughout Sec. Ill 


labels, though for al 
be regarded as fixed and arbitrary.) 


we do not suppress any of the mode 
intents and purposes l and m can 


An essential element in solving (3.1l is to obtain inde¬ 
pendent homogeneous solutions to (3.8), either through 
numerical integration (after setting causal boundary con¬ 
ditions at r* —i oo and r* -A — oo) or through use of an¬ 
alytic function (MST) expansions [6T]. We denote these 
unnormalized solutions as A^ n (r), where 

*Ln(r* + 00 ) - e™', 

Xr mn (n -A -oo) ~ e—*. 


(3.9) 

(3.10) 


A Green function is formed from these two linearly inde¬ 
pendent solutions and integrated over the source function 
Zimn(r) to obtain the particular solution of (|3.8[) 


Si m (t,r) = Gi m (t)S[r-r p {t)] 

+ F lm (t)S'[r-r p (t)]. [ ' j 

The time dependent functions Gi m (t) and Fi m (t) arise 
from tensor spherical harmonic decomposition m of the 
stress-energy tensor of the point mass and enforcement of 
the delta function constraints r -A r p {t) and ip -A ip p (t). 
Like the potential, their form depends upon parity. For 
l + Til even we use the Zerilli-Moncrief source, and for 
l + m odd we use the Cunningham-Price-Moncrief source 
(see m for details). 

As explained in Sec. |TlJ the eccentric motion of the 
source is characterized by two fundamental frequencies, 
and Q r . As such, we can represent the master func¬ 
tion and the source as Fourier series 


Ximn(r) = c+ mn (r) X+ mn (r) + c lmn {r) X lmn {r), 

where the normalization functions in the source 
are given by the integrals 

1 f r dr' - _ 

C Ln( r )=w - / J7ZF X lmn( r ') Z lmn(r'), 

vy lmn J r min J V ) 

1 f V max 

= Jr J(y) X ^nn( r ) Z l™n(r )■ 

Here Wi mn is the Wronskian 


(3.11) 

region 


(3.12) 


Wl mn 


f( r ) 



dX+ 


dr 


v+ 

yY lmn 



(3.13) 


OO 


4 f /m(Lr) = 

Xi mn (r)e~ iut , 

n=— oo 

(3.3) 

Sim(t,r) = 

oo 

Z lmn (r)e~ iut , 

(3.4) 


n =—oo 


with the mode frequencies being functions of both fun¬ 
damentals 


ui = u> mn = mflp + n£l r , m,n£l (3.5) 
The series coefficients are formally found by integrating 


While the expression in Eqn. ( |3. 11 ) is indeed a solution 
to Eqn. (3.8), it is not ideal. The singular nature of the 
TD source (3.2) results in Gibbs behavior in the Fourier 
synthesis (3.3) of '3q m at and near the particle location, 
leading to slow algebraic convergence. Exponential con¬ 
vergence can be restored by using the method of EHS, 
originally developed by Barack, Ori and Sago 1T51 . 

The first step in EHS is to extend the limits of integra¬ 
tion in (3.12) to include the full source region and obtain 
the normalization coefficients 


y;± _ 

V-y Imn 


l 

W^lmn 



^Imn( r )^™"( r ) 


/O') 


(3.14) 




























7 


These complex constants are in turn used to normalize 
the individual mode functions 


X lmn( r ) = G lmn X lmn 


0 )> 


(3.15) 


producing the FD EHS of Eqn. (3.8). Collectively, these 


normalized modes encode all the information about the 
source motion and are used to then define the TD EHS 




(f,r) = J2 X Ln( r ) e 


(3.16) 


As the FD EHS are each C °°, these Fourier sums con¬ 
verge exponentially for all r > 2 M. The sums are for¬ 
mally infinite in number, but in practice they are trun¬ 
cated once a specified accuracy is r eached. The desired 
particular TD solution to Eqn. ( |3.1| is then obtained by 
joining the outer and inner TD EHS: 

^lm(t,r) = ^tm e \ r - r pi t )] + ®Tm d l r p( t ) ~ r ] • ( 3 - 17 ) 

This weak solution can be computed everywhere, includ¬ 
ing the particle location, and it allows the metric and 
local gravitational self-force to be accurately determined 

M- 

There remains the practical issue of computing the 
G hnn■ For the RWZ problem, the source Zi mn (r) in 


Eqn. (3.14) is poorly behaved at the turning points be¬ 


cause of the presence of the S' term in (3.2) m- It was 
shown in that paper that the problem could be circum¬ 
vented by reversing the order of integration (see related 
examples in [351 ED])- To see this, s ubstit ute the Fourier 
transform integral for Zi mn (r) into (3.14) 


(7± _ 

w Zmn 


l 

Wi mn T r 


dr 


X Ln(r) 

f{r) 

x [ dt Sim{t, r)e lult 

J 0 


(3.18) 


Then subs titute for the TD source Si m (t,r) its singular 
form (3.2), exchange the order of integration, and inte¬ 


grate in r over the delta function terms. What remains 
of the calculation of Crf is an integral over time 


Cj =- 

^Imn 


i 


WimnTr 


+ 


r ^Inv 
JP 


,G tr 


2M_ jcT 1 dXf mn 


r 2 f 

p jp 


f P dr 


Fir 


(3.19) 


dt. 


The integrand is composed of obvious functions of time, 
such as Gim(t) and Fi rn (t). However, all of the other 
terms inside the square braces are now also functions 
of time, since the delta function maps r —> r p {t) [e.g., 
f P = f(r p (t)), X^ n (r) ->• Xf mn {r p {t))}. 

In summary, the RWZ BHP problem is solved by com¬ 
puting, for a sufficient range of l, m, and n, the inner 
and outer mode functions Xjf nn (r) (by ODE integration 
or analytic function expansion) and computing the inte¬ 


grals ( |3.19[ ) for the normalization coefficients Cff nn (using 


either IVP ODE integration m or a numerical quadra¬ 
ture routine 


B. SSI for the normalization coefficients 


SSI is a new modification in the way the normaliza¬ 
tion coefficients C Pmn are calculated. The key first step 
in developing SSI was actually the reversal in the order 
of integration described immediately above. The second 
essential step involves recognizing the periodic nature of 
the integrand in (3.19). The functions Fi m (t) and G; m (t), 

have complex time 


which contribute to the source Si n 
dependence because of the biperiodic motion and (typi¬ 
cally) incommensurate frequencies H r and The mo¬ 
tion in ip can be split into 


ip p (t) = VL^t + A 


(3.20) 


where the mean azimuthal advance is modulated by 
A ip(t), which is periodic in the radial motion. This <p p {t) 
enters source terms only through the spherical harmonic 
factor e ~ lm Vp{t ) ; w hi c h factors into: 

It is the mean azimuthal phase advance, at angular rate 
that makes source terms biperiodic. We can, how¬ 
ever, define functions Gi m and Fi m via 


Gim(t)=G lm (t)e imn 
F lm (t) = F lm {t) 


(3.21) 


that are strictly T r -periodic. Returning to Eqn. (3.19), 
we see that the factor, responsible for biperiod¬ 

icity, cancels with a corresponding factor from the Fourier 
transform kernel. We can replace the integral with 


Ct n = —[ r Etnit) e*' 
Wlmn-Lr J 0 




dt. 


where Ep^^t) are strictly T r -periodic functions 


Etnit) = T X L n Gl 


fp 


(3.22) 


(3.23) 


2 l dXl m \ 

„ , I -l In 


2f2 Imn f 
1 PJ P J 1 


dr 


The third, and most important, step toward SSI harks 
back to our earlier discussion in Sec. |IIC| of spectrally in¬ 
tegrating the orbit equations. There we showed that due 
to the C°° smoothness of (for example) dt p /dx = g{x) 
we could replace g(x) with an equally-spaced sampling 
gk — g(k Ay) of modest total number of samples N and 
achieve high-accuracy interpolation and integration. For 
source integration, the eq uivale nt step (to be justified 
momentarily) is to replace (3.22) with 


(7± _ 

w Zmn 


NWu 


N—l 

k—0 


In 


Mk 


inQrtk 


(3.24) 


where the time samples are tk = kT r /N , with k = 
0,... ,7V — 1. This remarkably simple sum is the heart 
of SSI. By replacing the integral in (3.22) with the sum 

























in (3.24), we avoid ODE integration and the calculation 


JV-l 


of the normalization coefficients is vastly sped up, open¬ 
ing the door to much higher accuracy applications [H]. 
What makes SSI work? Before we examine how well 


SSI performs, we first justify (3.24) as an appropriate re¬ 


placement for (3.22). The argument starts by noting the 
expec ted smoothness of the functions E Pmn (t) that enter 
(3.22). The contributing elements Fi m (t) and Gi m (t) are 
smooth C°° functions of the orbital motion. Similarly, 
the modes are smooth functions of r, and hence 

become smooth functions of time under the replacement 
r —¥ r p (t). Thus, for every Imn , the integrand in (3.22) 


is smooth and periodic. These properties suggest, just as 
they did in Sec. |IIB[ use of FS expansion. Indeed, the 


integral in Eqn. (3.22) looks like, under a cursory glance, 


the calculation of a set of FS coefficients. However, it is 
clear that C^ nn is not a spectrum of coefficients (in n) 
derived from a single function of time, but is instead cal¬ 
culated from a whole set (in n) of TD functions Ajfi n (f). 

Nevertheless, the Fourier series can be put to inves¬ 
tigative use and we introduce one for each Ififi n (t): 


Z7i=t 
^Imn 


(*>= E l , 


± 

Imnn' 


—in' fl r t 


n '=—oo 


with the coefficients given by 


1 


= -/ dtE* mn (t)e' 


i'Q r t 


(3.25) 


(3.26) 


If (3.25) is substituted in (3.22), and sum and integral 


are exchanged, we find that the normalization coefficients 


c ± = 
w /mn 


l 


Wl r 


"'Imnn 


(3.27) 


are proportional to the diagonal elements (n = n!) of the 
superset (over n and n') 
result is understandable: 


of FS coefficients fit, , 


The 

the integral in (3.22) simply 
picks out the nth harmonic in the nth function Ef^ (t). 

To complete the argument, we may assume (and nu¬ 
merically verify) that the smoothness of a source function 
EZ. (t) implies a rapidly falling (likely geometric) spec- 


lmn\ 

trum for £, 


Imnn' 


as n! —> ±oo. As we argued in Sec. 


IIB 


for any given accuracy goal, this implies the spectrum 
can be truncated at some sufficiently negative and posi¬ 
tive values of n'. Truncation, in turn, means that we have 
replaced the original source function with a bandlimited 
approximation. Bandlimiting then argues for replacing 
the source function (yet again), this time with a set of 
discrete, equally-spaced samples E^ nn (tk)- Because the 
source function is periodic, the discrete sampling is finite 
in number (say N). We can then use the DFT to relate 
the discrete sampling representation of the source to a 
discrete, finite spectrum (and vice versa) 


E Ln^k) = 


N—l 

V e ± 

/ j Imnn' 
n'—O 


— in' Q, r t k 


(3.28) 


£Lnn’ = x E E Ln(t*) ■ 


(3.29) 


k —0 


The DFT spectrum fifi nn / is distinct from the FS spec¬ 
trum £{? nnn / ; and the former will display periodicity in 
the FD, fi± , +jN = fi tun" However, for sufficiently 
large N and between the negative and positive Nyquist 
frequencies, the two spectra can be made nearly indistin¬ 


guishable. If we then set n! = n, replace fitmn in (3.27) 


with the DFT spectral component fifi ran , and s ubstitute 
into the same equation the DFT relation (3.29), we have 
derived our SSI formula Eqn. (3.24). 


We can provide a summary of this discussion, and the 
derivation, through a sequence of replacements: 


Cf = - 1 - f r dt Ef (t) e i: 

Imn tt t rji I Imn v / 

Wlmn-Lr J 0 




(3.30) 


1 


J 0 

1 


dte^t SLnn’ 


n' ——oo 


W^lmnTr J q 


- 1 — f 

Wlm n T r Jq 


1 


dte^ £ £t nn .e 


—in'Q r t 


dt e inQrt ^ c± 


E 


Imnn' ^ 


—in'£l r t 


W^lmn^r Jq 


dt e mne 


—in'Cl r t 


n—n. 
N—l 


x-YEf (t k ) e in ' Qrtk 

j\T / j lmn\ 


k—0 


N—l 


NWi, 


E E^E(4)e- 


Q r t k 


$nn' 


N—l 


NW, 


Imn 


E e 


inCl r t k 


k=0 


The two approximate (but typically spectrally accurate) 
steps are indicated. 

What is involved in practical use of SSI? Another way 
of asking this question is: if we make N discrete sam¬ 
ples of each source function and sum them in ( 3.24 1, for 
which and how many n’s should we compute C lmn l We 
do not presently have an exact answer, but we have an 
effective, practical procedure. To see the issue, consider 
Fig-0 There we show gravitational wave energy fluxes 
per harmonic n at r = oo for the l = 2, m = 2 mode 
(essentially proportional to |C^ n | 2 ). We might expect, 
for a given N, to begin near n = 0 and see a spectrum 
that descends on either side until hitting a Nyquist point 
(at about n = ±N/ 2). That is roughly, but not exactly, 
what is observed. The problem is that C^ nn is not, as 


a function of n, a DFT spectrum. If we consider (3.27), 
clearly the Wronskian Wi mn should not be expected to 






























9 




FIG. 4. Aliasing effect from oversampling SSI in the FD. 
Shown here are energy flux data from an orbit with p = 10 20 
and e = 0.01 for modes with l = 2, m = 2. The energy flux 
from successive n modes falls off exponentially when com¬ 
puted away from the peak harmonic. Note that one harmonic 
(n = —m = —2) is nearly static, which decreases its flux by 
more than 100 orders of magnitude. As higher positive and 
negative n are computed, the fluxes reach Nyquist-like notches 
and oversampling in n beyond those points leads to increases 
in flux similar to aliasing. The locations of the minima scale 
with but are not equal to ±A/2. 


display a periodicity in n. Even the DFT spectra, while 
having the periodicity in ri, ^ mn n , +jN = will 

not have a periodicity in the diagonal elements £^ nnn as 
a function of n. Nevertheless, if we sample in n 

for \n\ > N/2 we observe a succession of Nyquist-like 
notches and peaks, similar to aliasing in the DFT but 
without exact periodicity. From a practical standpoint, 
we compute and use the spectrum in n down to the first 
Nyquist-like notch on each side and calculate no further. 
The code marches forward on each side, finds the minima, 
and discards contributions beyond those points. 

Fig. a shows this aliasing phenomenon. There we de¬ 
liberately compute and display energy fluxes for a few 
harmonics beyond the first Nyquist notch on each side 
of the central maximum. We show the same fluxes com¬ 
puted with four different spectral resolutions. The ex¬ 
ponential fall in the spectrum is evident. These calcula¬ 
tions were made possible not only by use of SSI but also 
Mathematica’s arbitrary precision arithmetic. As N be¬ 
comes larger, we approach the FS, or continuum, limit. 
It is clear from the vantage point of high resolution that 
the best thing to do at lower resolution is halt the mode 
calculations at the Nyquist notches. This assumption is 
borne out by considering Fig. [5j This figure displays the 
differences in fluxes between those computed at resolu¬ 
tions of N = 40, 60, 80 and those found with N = 100. 
The error in the discrete representation is well bounded 
in the region between the first Nyquist points by the max¬ 
imum error at one of the notches. 

Operating at high accuracy (e.g., 200 decimal places), 


FIG. 5. SSI at high accuracy. Absolute differences (errors) 
are shown in self-convergence tests. The same data found 
in Fig. [4] are used to compute differences (per harmonic) in 
the fluxes between the three lower resolutions and the highest 
(N = 100) resolution. For each of the three lower resolutions 
(N = 40, 60, 80) the errors are well bounded by the accuracy 
criteria set by errors at the Nyquist notches. 


the MST code makes a prediction of how large N needs 
to be in order that the Nyquist notches lie (just) below 
the specified error level. We have observed that adequate 
sampling for SSI is always sufficient for comparably ac¬ 
curate orbit integration. The prediction for N is tested 
and if mode fluxes do not reach the error level at the 
Nyquist notches, then a new value of N is chosen and 
the calculation is repeated. In the MST code, the now 
vastly reduced number of function evaluations can still 
be expensive if N is set too generously. It is important 
to note th at the key formula for SSI, Eqn. (3.24) (or more 
properly (3.31)—see below), is an 0(N 2 ) procedure, and 
so it is essential to find a near minimum value of N for 
a given accuracy goal. 

How well does SSI work? We have demonstrated nu¬ 
merically in Figs. [2] and [5] the presence of exponential 
convergence. It is not that the gravitational wave fluxes 
fall geometrically (a known result), but that the gap be¬ 
tween resolutions (i.e., error in substituting the DFT for 
the Fourier series) falls exponentially with increases in N. 
We can, however, go a step further and make the rate of 
exponential convergence even faster by introducing one 
final modification. 

SSI is exponentially convergent because the periodic 
functions, E^ nn {£), being sampled are C°°. However, 
there is no requirement that the periodic motion be de¬ 
scribed by t. Any C°° reparametrization t —> A (t) should 
be expected to also give rapidly convergent sums. This is 
true, for example, in switching from t to the relativistic 
anomaly x- We have found empirically, though, that use 
of x as the curve parameter substantially improves the 
rate of exponential convergence. 

To effect this change we rewrite Eqn. (3.22) with x as 
the independent variable. Then the periodic motion is 
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FIG. 6. Comparison between SSI with y sampling and t sam¬ 
pling as a function of N and eccentricity. The two sampling 
schemes are tested by examining convergence of the l = 2, 
m = 2, n = 0 energy flux using the MST code. All orbits 
have p = 10 3 . 


divided into equally spaced steps Ay = 27 t/N, the inte¬ 
grand is discretely sampled, and the integral is replaced 
with the sum 


(j± 
w Imn 


NW lmn 


dt n - 


Y,dx E Lnlt(x k )]e inn ^. (3.31) 


As Fig. [6] demonstrates, substantially fewer y-samples 
are required than f-samples for SSI to reach a prescribed 
error level. This is especially true of high eccentricity 
orbits. As a practical matter, it is also easier to find ip p 
and r p evenly sampled in % than in t. Finally, we com¬ 
ment that it was merely a hunch (though one informed 
by experience with the problem) that y might provide 
a better measure and more rapid convergence. It is an 
open question whether there is another parametrization 
of the orbit that yields even faster convergence. 

The numerical SSI results shown so far have involved 
the high accuracy MST code. But SSI also aids in dou¬ 
ble precision C coded calculations. Its benefit is shown 
clearly in Fig. [TJ where we mark the accuracy reached in 
computing a normalization coefficient (C^ 0 ) as a func¬ 
tion of the number of source term evaluations, which 
serves as a proxy for computational load. We compare 
SSI to an IVP ODE integration using an 8th order Runge- 
Kutta routine. SSI has exponentially converging accu¬ 
racy with increases in function calls (i.e., increases in N). 
In contrast, the Runge-Kutta routine, with its algebraic 
convergence, struggles to reach high accuracies. 

Do we really need SSI? The answer depends upon the 
application. At double precision the answer is clearly no, 
but SSI is likely much more efficient (and hence faster). 
The real critical requirement for SSI comes in high accu¬ 
racy eccentric orbit calculations. Consider Fig. [7] again 
and the scaling of ODE integration. At an accuracy goal 
of 200 digits, even an efficient algorithm like 8th order 


Runge-Kutta would take of order 10 22 steps to integrate 
through an eccentric orbit source region! Without SSI or 
a comparable spectral method, these calculations simply 
cannot be done. 


C. SSI and the midpoint and trapezoidal rules 


We developed SSI with the convergence of Fourier se¬ 
ries and concepts in signal processing (e.g., sampling the¬ 
orem, use of the DFT/FFT, etc) firmly in mind. The ap¬ 
plication of SSI to orbit integration does in fact simply 
use the DCT, a special case of the DFT. In source inte¬ 
gration, even though the key formula, ( |3.31 ) or (3.24|, is 
not a DFT, we used the DFT to provide an understand¬ 
ing of the rapid convergence of the sum. The essential 
point was to see that rapid convergence in the FD with 
a modest number N of spectral elements could translate 
into representing the behavior in the TD with equally 
modest sampling. If the representation has sufficient ac¬ 
curacy, then interpolation and integration can be made 
accurately as well. 

Yet, if we step back and examine the sums [(3.31) or 
(3.24)] that we use in SSI, a curious fact jumps out: 
they appear to be nothing more than simple Riemann 
sums. Given the sampling, their form appears to be a 
use of the left rectangle rule. However, with the inher¬ 
ent periodicity in y, the left rectangle rule is equivalent 
to the trapezoidal rule and, with a half interval shift in 
the equal-sized y bins, it is also just the midpoint rule. 
But these are just the lowest-order approximations for 
an integral, with error bounds, 0(1/N 2 ), that are alge¬ 
braic in the number of divisions of the interval! How can 
their use be giving a vastly faster rate of convergence? 
The answer lies in the periodicity and smoothness of the 
summands. After developing the method we came upon 
a paper [BT l that discusses this surprising behavior in 
other contexts and nicely provides a set of example cal¬ 
culations. A more recent and exhaustive discussion is 
found in [55]. In black hole perturbation work, Fujita, 
Hikida, and Tagoshi [52] made use of the trapezoidal 
rule for source integration, but did not explicitly note 
or demonstrate the exponential convergence or push the 
results beyond double precision. The remarkably rapid 
convergence of the trapezoidal rule in special cases is ap¬ 
parently well known in certain numerical analysis circles, 
though it also appears to be something that is constantly 
being rediscovered ever since Poisson’s original finding in 

1827 55]. 


Since we have shifted the viewpoint momentarily to 
thinking about Riemann sums and quadrature formu¬ 
lae, what about higher-order methods like Simpson’s 
rule? Since Simpson’s rule generally has a stronger error 
bound (implying presumably faster convergence) than 
the trapezoidal rule, might its use in SSI allow us to con¬ 
verge even faster? Alas, the answer is no, as a quick test 
demonstrated. A discussion and example can be found 
in [ 54] . 
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FIG. 7. Efficiency of spectral source integration in comparison to ODE integration in a RWZ application at double precision. 
The RWZ normalization constant C22 0 is computed for various eccentricities in orbits with p = 10. The ODE integration uses 
the Runge-Kutta-Prince-Dormand 7(8) [62] routine rk8pd of the GNU Scientific Library (GSL) 1631 . 


IV. SPECTRAL SOURCE INTEGRATION IN 
LORENZ GAUGE 

After developing SSI for use in the RWZ formalism, 
we turned our attention to Lorenz gauge and were able 
to successfully apply the method to coupled systems of 
equations. Lorenz gauge breaks down into several sys¬ 
tems of different orders that depend on (1) parity, (2) 
mode order (either low l = 0,1 or high l > 2), and (3) the 
special static (u> = 0) case. See [29] and !30| for details. 
Here we simply demonstrate the principles of incorporat¬ 
ing SSI by focusing only on the odd-parity equations. 


A. Odd-parity Lorenz gauge and EHS 


In Lorenz gauge, odd-parity perturbations can be de¬ 
scribed by the two amplitudes, h l t m and hf™, with the 
third, h l ™ , obtainable from the gauge condition. The 
reduced-order coupled system in the TD is 




d,h l : 


r* r 


(4.1) 


2p + (l + 2)(l-l)f . 


K m = f 2 PL, 


□ hi m + 2 ^ r drhl m - 2(r n 3M) dth l r (4.2) 


2 ^r'^r 

r-z 

c2 


r* r 2 f 


2p + (l + 2)(l-l)f, 


hT = —P, 


Im • 


It is convenient in what follows to write the fields and 
their sources in a vector notation 


The Lorenz gauge source terms are proportional to the 
delta function 5[r — r p (t )] (in contrast to RWZ gauge 
where the source also has a S' term), allowing the source 
vector to be expressed in terms of a time-dependent vec¬ 
tor amplitude Vi m (t), 


Vi m (t,r) = Vi m (t)S[r - r p (t)]. (4.4) 

The field Bi m and source V; m can be e xpre ssed as 
Fourier series, analogous to Eqns. (3.3) and (3.4), 


OO 


Bi m {t,r) = 

£ B lmn (r)e- iut , 

71— — 00 

(4.5) 

Vim(t,r) = 

•1 

T 

s 

s 

8 w 

(4.6) 


n =—00 


The FS coefficients are formally found via the integrals 


&imn{r) = — / dt Bi m (t,r)e lult , (4.7) 


L T J 0 


V Imn{? ) — 


T, 


dt Vim{t,r) e u 


(4.8) 


r JO 


With these definitions, we henceforth in this section sup¬ 
press the TD mode labels Im and the FD labels Imn 
whenever no confusion might arise. However, for clarity 
we attach a tilde to denote FD quantities. 

In the FD, the field equations (4.1) and (4.2) take the 
following form 


Bim(t, r) 


h[ m 

fh l r 


Vi m (t,r) 


PPL 

-fPL 


(4.3) 


d^B + Cd r ,B + T>B = V. (4.9) 


The matrices C and D that couple the equations are 
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given by 


c = i 


-M 0 
0 r - 3 M 


v = [ ^_vi±l±m^r n 


(4.10) 


2 iui 
„2 


0 -M 
r — 3M 0 


where I is the 2x2 identity matrix. 

The EHS method carries over to Lorenz gauge EH EDI. 
Four fundamental independent homogeneous solutions to 


~± 


Eqn. (4.9) are denoted by B { , with i = 0,1. The ± su¬ 


perscript delineates causal asymptotic behavior, with + 
indicating an outgoing wave at r* = oo and — indicating 
a downgoing wave at r* = — oo. A Green function con¬ 
structed from these arbitrarily normalized modes yields 
the solution to the inhomogeneous system 


B = Bq Cq + B[ cf + Bq Cq + By Cy , 


(4.11) 


once the c~(r) are determined by integrating the first- 
order linear system 


M(r) 


clrv Cn 


Cq" 

d Ttt Cy J 


0 

V(r) 


(4.12) 


Here M is the (Imn dependent) 4x4 Wronskian matrix 


find the normalization coefficients with an integral over 
time, 


—Cn 


-Cy 

Cq 

Ct 


1 

f Tr 1 1 

0 


/ —Mr 1 


Tr 

o 

v(t) _ 


e iuit dt. 


(4.16) 


In this last equation the script p indicates time depen¬ 
dence via the mapping r r p (t). With the coefficients 
available, the FD and TD EHS (respectively) are con¬ 
structed 


6 ( r ) = CqB 0 (r) + CfBy (r), 


B ± (t,r)= ^(r) 


(4.17) 

(4.18) 


The solution to the system in the TD, Eqns. (4.1) and 
(4.2), is then 


B(t, r ) = B + 9 [r — r p (t)] + B 9 [ r p (t) — r\ . (4.19) 

T he ke y to EHS in Lorenz gauge is solving systems 
like (4.16) for the normalization coefficients. In previous 
work WED] these equations were treated as IVPs and 
solved with ODE integration. That numerical approach 
can be replaced with SSI to achieve spectral convergence, 
as we outline next. 


B. Spectral source integration for odd-parity 
normalization constants 


M(r) 


B 0 By Bq By 

d r ,Bo d r ,By dr,Bo dr,By \ ’ 


(4.13) 


and 0 is the rank = 2 column vector. 

Solving for the functions cf(r) can be avoided through 
use of the method of EHS. Instead, as in Sec. m we 
solve for normalization coefficients, construct FD EHS 
and then TD EHS, and thus circumvent producing Gibbs 
behavior in the source region and at the particle location. 
For the system at hand, we define the normalization con¬ 
stants Cf as 


The Lorenz gauge employment of SSI is virtually iden¬ 
tical to RWZ gauge. As in Eqn. (3.21), we can extract 
from the biperiodic source term v(t) the piece that is 
periodic in T r by defin ing v(t) = v(t)e lmn ‘ f,t . Once sub¬ 
stituted in Eqn. (4.16) we find 


-Co 

-Cy 

c£ 

Ct 


1 

Tr 


E{t) e 


inQ r t 


dt , 


(4.20) 


Ct=4(r n 


<)) C i — Cj (r m ; n ), 


(4.14) where we define the vector 


and obtain them via the integrals 

'-Co 

-Cy 

Ct 

. ct 

In the expression above, M(r) -1 is the inverse of the 
Wronskian matrix M(r). We next insert the integral ex¬ 
pression for V(r), reverse the order of integration, and 


f"T max 1 

0 

L. /M Mlr:i 

,V(r)_ 


dr. (4.15) 


m = jMp 1 

Jp 


0 

v(t) 


(4.21) 


Both E(t) and e m! * rt are periodic in T r . The vector E(t), 
which depends on the FD labels Imn , is the equivalent 

of Etmnit in E q n - d 3 23 D - 

The logical steps in implementing SSI carry over from 
Sec. lfflBl 


• The vector E(t) (carrying labels Imn ) consists of 
periodic, C°° functions. 
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• Each can be represented as a Fourier series with 
spectrum £ n i with n' —> ±oo. 

• The Fourier series spectrum can be truncated to 
some n' min < n! < n' mSLK subject to an accuracy 
goal. 

• The approximate but very accurate truncated 
Fourier series is a bandlimited function. 

• The Nyquist-Shannon sampling theorem implies 
the truncated Fourier series representation can it¬ 
self be replaced in the TD with discrete sampling. 

• Sampling plus periodicity implies a discrete repre¬ 
sentation of finite length N. 

• Finite sampling representation in the TD implies 
one-to-one correspondence via the DFT with a FD 
periodic spectrum £„/. 

• The DFT spectrum within the first Nyquist minima 
approximates well the original Fourier series spec¬ 
trum if N is sufficiently large, allowing £ n i —> £ n >. 


Based upon this chain of reasoning, the integral (4.20) 
can be replaced with an exponentially convergent sum 


-Co 

-c'r 

Co + 

. c+ 


1 

N 


N-l 


Y E(t k ) e mnrtk 

k=0 


(4.22) 


where again t k = kT r /N. This is SSI for the systems of 
equations found in Lorenz gauge. 

Even parity involves a larger linear system. The matrix 
inversion in evaluating E(t k ) at the sample points is the 
most expensive task in double precision application. We 
use LU decomposition and take advantage of the symme¬ 
try M(tfc) = M (T r — t k ), so that LU decompositions of 
M are only necessary at N /2 + 1 points. 

As in Sec. |III B[ we also attain a higher rate of expo¬ 
nential convergence by switching from t parametrization 
to x■ The adjustment to Eqn. (4.22) is straightforward 


-Co 

-cr 

C 0 + 

. c+ 


N 


N-l 


E 


dt p 

dx 


E[t(Xk )] 


inO. T t(xk) 


(4.23) 


where as before Xfc = ^irk/N. 

We have implemented SSI in just this way as a modi¬ 
fication of the Lorenz gauge code described in Ref. [50] . 
SSI is particularly beneficial in Lorenz gauge, where a 
matrix must be inverted at each step (i.e., each function 
evaluation) in an integration. It is also beneficial that 
we know precisely where the sample locations are in the 


source region before computing the inner and outer ho¬ 
mogeneous solutions. The previous method found the 
normalization coefficients by integrating a large simulta¬ 
neous system of ODEs through the source region (for the 
even-parity field this tallied to integrating 144 variables 
simultaneously). With prior knowledge of the sample lo¬ 
cations, integration of the homogeneous solutions is de¬ 
coupled from the SSI for the normalization coefficients. 

As we described in Sec. |III B[ the number of sample 
points N is determined, iteratively if necessary, based on 
an error criterion. In all cases we have experience with, 
both RWZ and Lorenz gauge, it is the source integration 
(with SSI), not the spectral integration of the orbit, that 
sets the condition on N. In Lorenz gauge, the require¬ 
ment on N to meet a double precision error criterion in 
SSI is about a factor of 8 larger than required for a com¬ 
parably precise orbit integration. (With the MST code 
at a high accuracy of 200 digits, SSI requires an N that is 
about a factor of 2 larger than that required for compa¬ 
rably accurate orbit determination.) Because SSI shrinks 
so markedly the computational work in finding the nor¬ 
malization coefficients, our Lorenz gauge GSF code is 
sped up-overall by a factor of about 3 for eccentricities 
of order e ~ 0.7. 


V. CONCLUSIONS 

We have described in this paper a new method for 
achieving spectral accuracy and computational efficiency 
in calculating a broad class of black hole perturbation and 
gravitational self-force problems that entail generic or¬ 
bital motion. This class should include most problems in¬ 
volving a point-particle description of the small compact 
object and use of the background geodesics in a frequency 
domain calculation (i.e., geodesic self-force calculations). 
We have shown it applied both to the RWZ formalism (for 
individual master equations) and to Lorenz gauge (with 
coupled systems of equations) for eccentric binaries with 
a Schwarzschild primary. The method should extend to 
extreme-mass-ratio inspirals on Kerr as well, which will 
be addressed in subsequent work. Called spectral source 
integration (SSI), this method provides an exponentially- 
convergent calculation of the mode normalization coeffi¬ 
cients by replacing the more typically used ODE integra¬ 
tions in the source region. A simple modification of the 
underlying idea is also used to integrate the equations of 
orbital motion, to provide a consistent level of accuracy 
in determining source functions in the libration region. 

Use of SSI in double precision calculations will improve 
code speed and help ensure optimal accuracy. In con¬ 
trast, SSI is the sine qua non for calculating eccentric 
binaries using (MST) analytic function expansions at ex¬ 
traordinarily high accuracies (e.g., 200 decimal places). 
No algebraically convergent ODE solver will be able to 
calculate eccentric-orbit perturbations at hundreds of 
decimal places of accuracy. Any alternative to SSI will 
almost certainly be a similar technique using some other 
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spectral basis. A subsequent paper will describe use of 
SSI in an MST code to uncover new terms in the post- 
Newtonian expansion for eccentric binaries well beyond 
known 3PN order 35] , 


We then focus on the leading Newtonian term and seek 
to find its Fourier spectrum. Our derivation is similar to 
one found in ^4j. Adopting the notation 5 at(x) for the 
term in question, we first introduce complex exponentials 
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Appendix A: Exact Fourier spectrum for dt p /d\ as 
p —> oo 

Several figures in this paper have shown numerical ev¬ 
idence of exponential fall-off in the FD spectra of source 
functions in the perturbation and orbit equations. Here 
we demonstrate an exact calculation of the Fourier spec¬ 
trum in one simplified case. Consider the source function 


g(x) i 11 Eqn. (2.5) and make a post-Newtonian expansion 


dtp _ Mp 3 / 2 

d\~ 9X (1 + ecosy) 2 


The denominator in this expression can be factored 

_ 4Mp 3 / 2 <7 2 

9N{X ~ e 2 (l + ae ix ) 2 (! + ae~ ix ) 2 ’ 


(A3) 


by introducing 


1 


e (l - \A - e 2 ) , (A4) 

which is one of the roots of the quadratic equation 
a 2 — (2/e)<r + 1 = 0. We then make a partial fractions 
decomposition of Eqn. ([A3 1 


9n(x) = - 


Mp 3 / 2 


4 (l + cr 2 ) a 2 + 


2\„2 , 4ct 2 (<7 2 -1) 


e 2 (a 2 - l) 3 
8 ct 4 4ct 2 6 (cr 2 — l) 


(1 + ae lx ) 2 
8 cr 4 


1 + ae lx (l -f ae~ ix ) 1 + cre _IX 


(A5) 


Since it can be shown that |cr| < 1 for bound motion, each 
of the terms in Eqn. (A5| can be expanded in binomial 


or geometric series. The result is a Fourier series in y. 
Because of the symmetry of <7jv(x), the expansion reduces 
to a cosine series 


9n(x) = ^Go + ^Gn cos (ny), 

n—1 

and we find that the spectrum has the form 
8Mp 3 / 2 (-l) n+1 


(A6) 


Gn = 


e 2 (l — cr 2 ) 3 


[(n — l)cr n+4 — (n + l)cr" +2 ] . 


+ O (p 1 / 2 ^ . (Al) The exponential convergence of the series is evident. 


(A7) 


[1] A. Le Tiec, International Journal of Modern Physics D 
23, 1430022 (2014), arXiv:1408.5505 [gr-qc] 

[2] T. W. Baumgarte and S. L. Shapiro, Numerical Relativ¬ 
ity: Solving Einstein’s Equations on the Computer (Cam¬ 
bridge University Press, 2010). 

[3] L. Lehner and F. Pretorius, Annu. Rev. Astron. Astro- 
phys. 52, 661 (2014), arXiv: 1405.4840 [astro-ph.HE] 

[4] C. M. Will, Proceedings of the National Academy of Sci¬ 
ence 108, 5938 (2011), arXiv:1102.5192 [gr-qc], 

[5] L. Blanchet, Living Reviews in Relativity 17, 2 (2014) 
arXiv:1310.1528 [gr-qc] 

[6] S. Drasco and S. A. Hughes, Phys. Rev. D 69, 044015 

(2004) 


[7] L. Barack, Class. Quant. Grav. 26, 213001 (2009) 
arXiv:0908.1664 [gr-qc] 

[8] E. Poisson, A. Pound, and I. Vega, Living Rev. Rel. 14, 
7 (2011), arXiv:gr-qc/1102.0529 

[9] J. Thornburg, GW Notes, Vol. 5, p. 3-53 5, 3 (2011), 
arXiv: 1102.2857 [gr-qc] 

[10] A. Buonanno and T. Damour, Phys. Rev. D 59, 084006 
(1999), gr-qc/9811091 

[11] A. Buonanno, Y. Pan, H. P. Pfeiffer, M. A. Scheel, L. T. 
Buchman, and L. E. Kidder, Phys. Rev. D 79, 124028 
(2009), arXiv:0902.0790 [gr-qc] 

[12] T. Damour, Phys. Rev. D 81, 024017 (2010) 

arXiv:0910.5533 [gr-qc] 






















15 


[13] T. Hinderer, A. Buonanno, A. H. Mroue, D. A. Hem- 
berger, G. Lovelace, H. P. Pfeiffer, L. E. Kidder, M. A. 
Scheel, B. Szilagyi, N. W. Taylor, and S. A. Teukolsky, 
Phys. Rev. D 88, 084005 (2013), arXiv: 1309.0544 [gr-qc] 

[14] T. Darnour, ArXiv e-prints (2013), arXiv:1312.3505 [gr- 
qc] 

[15] A. Taracchini, A. Buonanno, Y. Pan, T. Hinderer, 

M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace, 
A. H. Mroue, H. P. Pfeiffer, M. A. Scheel, B. Szilagyi, 

N. W. Taylor, and A. Zenginoglu, Phys. Rev. D 89, 
061502 (2014), arXiv: 1311.2544 [gr-qc] 

[16] Y. Mino, M. Sasaki, and T. Tanaka, Phys. Rev. D 55, 
3457 (1997). 

[17] T. C. Quinn and R. M. Wald, Phys. Rev. D 56, 3381 
(1997), 

[18] P. Amaro-Seoane, J. R. Gair, A. Pound, S. A. 
Hughes, and C. F. Sopuerta, ArXiv e-prints (2014), 
arXiv: 1410.0958 

[19] D. A. Brown, J. Brink, H. Fang, J. R. Gair, C. Li, 
G. Lovelace, I. Mandel, and K. S. Thorne, Physical Re¬ 
view Letters 99, 201102 (2007), gr-qc/0612060 

[20] P. Amaro-Seoane, J. R. Gair, M. Freitag, M. C. Miller, 
I. Mandel, C. ,J. Cutler, and S. Babak, Class. Quant. 
Grav. 24, R113 (2007), arXiv:astro-ph/0703495 

[21] L. Barack and A. Ori, Phys. Rev. D 61, 061502 (2000) 
arXiv:gr-qc/9912010 

[22] K. Martel, Phys. Rev. D 69, 044025 (2004) 

[23] L. Barack and N. Sago, Phys. Rev. D 75, 064021 (2007) 
arXiv:gr-qc/0701069 

[24] S. E. Field, J. S. Hesthaven, and S. R. Lau, Class. Quant. 
Grav. 26, 165010 (2009), arXiv:0902.1287 [gr-qc] 

[25] P. Canizares and C. F. Sopuerta, ArXiv e-prints (2014), 
arXiv:1406.7154 [gr-qc] 

[26] S. Detweiler, Phys. Rev. D 77, 124026 (2008) 

arXiv:0804.3529 [gr-qc] 

[27] S. Hopper and C. R. Evans, Phys. Rev. D 82, 084010 

( 2010 ) 

[28] S. Akcay, Phys. Rev. D 83, 124026 (2011) 

[29] S. Akcay, N. Warburton, and L. Barack, Phys. Rev. D 
88, 104009 (2013), arXiv: 1308.5223 [gr-qc] 

[30] T. Osburn, E. Forseth, C. R. Evans, and S. Hopper, 
Phys. Rev. D 90, 104031 (2014) 

[31] I. Vega and S. Detweiler, Phys. Rev. D 77, 084008 (2008) 
arXiv:0712.4405 [gr-qc] 

[32] I. Vega, B. Wardell, and P. Diener, Classical and Quan¬ 
tum Gravity 28, 134010 (2011), arXiv:1101.2925 [gr-qc] 

[33] B. Wardell, I. Vega, J. Thornburg, and P. Diener, Phys. 
Rev. D 85, 104044 (2012), arXiv:1112.6355 [gr-qc] 

[34] M. Casals, S. Dolan, A. C. Ottewill, and B. Wardell, 
Phys. Rev. D 79, 124043 (2009), arXiv:0903.0395 [gr-qc] 

[35] M. Casals, S. Dolan, A. C. Ottewill, and B. Wardell, 
Phys. Rev. D 88, 044022 (2013), arXiv: 1306.0884 [gr-qc] 

[36] B. Wardell, C. R. Galley, A. Zenginoglu, M. Casals, S. R. 
Dolan, and A. C. Ottewill, Phys. Rev. D 89, 084021 
(2014), arXiv:1401.1506 [gr-qc] 

[37] N. Sago, L. Barack, and S. L. Detweiler, Phys. Rev. D 
78, 124024 (2008), arXiv:0810.2530 [gr-qc] 

[38] L. Blanchet, S. Detweiler, A. Le Tiec, and B. F. Whiting, 
Phys. Rev. D 81, 064004 (2010), arXiv:0910.0207 [gr-qc] 

[39] L. Blanchet, S. Detweiler, A. Le Tiec, and B. F. Whiting, 
Phys. Rev. D 81, 084033 (2010), arXiv: 1002.0726 [gr-qc] 

[40] R. Fujita, Progress of Theoretical Physics 128, 971 
(2012), arXiv:1211.5535 [gr-qc] 


[41] A. G. Shah, J. L. Friedman, and B. F. Whiting, Phys. 
Rev. D 89, 064042 (2014), arXiv: 1312.1952 [gr-qc] 

[42] A. G. Shah, Phys. Rev. D 90, 044025 (2014) 
arXiv: 1403.2697 [gr-qc] 

[43] A. G. Shah and A. Pound, (2015), arXiv:1503.02414 [gr- 
qc] 

[44] S. Mano, H. Suzuki, and E. Takasugi, Progress of The¬ 
oretical Physics 95, 1079 (1996), gr-qc/9603020 

[45] E. Forseth, C. R. Evans, and S. Hopper, (in preparation, 
2015). 

[46] C. Hopman and T. Alexander, Astrophysical Journal 
629, 362 (2005), arXiv:astro-ph/0503672 

[47] B. Moore, M. Favata, K. Arun, and C. Mishra (2015) , 
Talk presented at the APS April meeting held in Balti¬ 
more, MD, http://meetings.aps.org/Meeting/APR15/ 
Session/R13.1 

[48] C. Cutler, D. Kennefick, and E. Poisson, Phys. Rev. D 
50, 3816 (1994) 

[49] L. Barack, A. Ori, and N. Sago, Phys. Rev. D 78, 084021 
(2008), arXiv:0808.2315 

[50] N. Warburton and L. Barack, Phys. Rev. D 83, 124038 
(2011), arXiv:1103.0287 [gr-qc] 

[51] S. Drasco and S. A. Hughes, Phys. Rev. D 73, 024027 
(2006), gr-qc/0509101 

[52] R. Fujita, W. Hikida, and H. Tagoshi, Prog. Theor. Phys. 
121, 843 (2009), arXiv:0904.3810 [gr-qc] 

[53] L. Barack and N. Sago, Phys. Rev. D 81, 084021 (2010), 
arXiv: 1002.2386 [gr-qc] 

[54] C. Darwin, Proceedings of the Royal Society of London. 
Series A: Mathematical and Physical Sciences 263, 39 
(1961). 

[55] I. S. Gradshteyn, I. M. Ryzhik, A. Jeffrey, and D. Zwill- 
inger, Table of Integrals, Series, and Products, Seventh 
Edition by I. S. Gradshteyn, I. M. Ryzhik, Alan Jef¬ 
frey, and Daniel Zwillinger. Elsevier Academic Press, 
2007. ISBN 012-373637-4 (2007). 

[56] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and 
B. P. Flannery, Numerical Recipes in C : The Art of Sci¬ 
entific Computing, 2nd ed. (Cambridge University Press, 
Cambridge, UK, 1993). 

[57] R. Fujita and W. Hikida, Classical and Quantum Gravity 
26, 135002 (2009), arXiv:0906.1420 [gr-qc] 

[58] F. Stenger, SIAM Review 23, pp. 165 (1981) 

[59] N. Ahmed, T. Natarajan, and K. R. Rao, IEEE Trans. 
Comput. 23, 90 (1974) 

[60] K. R. Rao and P. Yip, Discrete Cosine Transform: Algo¬ 
rithms, Advantages, Applications (Academic Press Pro¬ 
fessional, Inc., San Diego, CA, USA, 1990). 

[61] M. Sasaki and H. Tagoshi, Living Reviews in Relativity 
6, 6 (2003), gr-qc/0306120 

[62] P. Prince and J. Dormand, Journal of Computational and 
Applied Mathematics 7, 67 (1981) 

[63] “Gnu scientific library,” http: //www. gnu. org/ 
software/gsl/ 

[64] J. A. C. Weideman, American Mathematical Monthly 
109, 21 (2002). 

[65] L. N. Trefethen and J. Weideman, SIAM Review 56, 385 
(2014). 

[66] S. D. Poisson, Memoires de l’Academie Royale 
des Sciences de PInstitut de France 4, 571 (1827), 

http://hdl.handle.net/2027/mdp.39015077785494? 
urlappend="/,3Bseq=769 





